<!DOCTYPE html>
<html lang="en">
<head>
        <title>Theremin Simulation</title>
        <meta charset="utf-8" />
	<meta name="viewport" content="width=device-width, initial-scale=1.0" />
        <link rel="shortcut icon" href="../theme/images/favicon.ico"/>
        <link rel="stylesheet" href="../theme/css/main.css" type="text/css" />
        <link href="www.j-marjanovic.io/feeds/all.atom.xml" type="application/atom+xml" rel="alternate" title="j-marjanovic.io Atom Feed" />

        <!--[if IE]>
                <script src="http://html5shiv.googlecode.com/svn/trunk/html5.js"></script><![endif]-->

        <!--[if lte IE 7]>
                <link rel="stylesheet" type="text/css" media="all" href="../css/ie.css"/>
                <script src="../js/IE8.js" type="text/javascript"></script><![endif]-->

        <!--[if lt IE 7]>
                <link rel="stylesheet" type="text/css" media="all" href="../css/ie6.css"/><![endif]-->
<script src="https://ajax.googleapis.com/ajax/libs/jquery/1.8/jquery.min.js" type="text/javascript"></script>
<script type="text/javascript" src="theme/jquery.cookiesdirective.js"></script>

</head>

<body id="index" class="home">
<script type="text/javascript">
	// Using $(document).ready never hurts
	$(document).ready(function(){

		// Cookie setting script wrapper
		var cookieScripts = function () {
			// Internal javascript called
			console.log("Running");
		
			// Loading external javascript file
			$.cookiesDirective.loadScript({
				uri:'external.js',
				appendTo: 'eantics'
			});
		}
	
		/* Call cookiesDirective, overriding any default params
		
			*** These are the defaults ***
				explicitConsent: true,
				position: 'top',
				duration: 10,
				limit: 0,
				message: null,				
				cookieScripts: null,
				privacyPolicyUri: 'privacy.html',
				scriptWrapper: function(){},	
				fontFamily: 'helvetica',
				fontColor: '#FFFFFF',
				fontSize: '13px',
				backgroundColor: '#000000',
				backgroundOpacity: '80',
				linkColor: '#CA0000'
				
		*/
		
		$.cookiesDirective({
			privacyPolicyUri: 'myprivacypolicy.html',
			explicitConsent: false,
			position : 'bottom',
			scriptWrapper: cookieScripts, 
			cookieScripts: 'Google Analytics', 
			backgroundColor: '#52B54A',
			linkColor: '#ffffff'
		});
	});
</script>

	
  <!-- <header id="banner" class="body"> -->
  <!--               <h1><a href="../"><img src="http://www.launchyard.com/images/logo.png" /></a></h1> -->
  <!--       </header> --> 

  <div class="LaunchyardDetail" style="align:right;">
    <!-- <p> -->
    <!-- <img src="../theme/images/blue-pin.png" width="100" height="100" alt="Graph icon"> -->
    <!-- </p> -->
    <p><a id="sitesubtitle" href="../">j-marjanovic.io</a></p>

	<br>
    <p style="float: right; margin-right: 50px;"><a id="aboutlink" href="../pages/about.html">About</a></p>

    <br>
	<p style="float: right; margin-right: 50px;"><img src="../theme/images/icons/rss.png"> <a id="aboutlink" href="../feeds/jan-marjanovic.atom.xml">Atom feed</a></p>
    <br>

  </div>

<section id="content" >
    <div class="body">
      <article>
        <header>
          <h1 class="entry-title">
            <a href="../drafts/theremin-simulation.html" rel="bookmark"
               title="Permalink to Theremin Simulation">Theremin Simulation</a></h1>

        </header>

        <div class="entry-content">
<div class="post-info">
	<ul>
        <li class="vcard author">
                 by&nbsp;<a class="url fn" href="../author/jan-marjanovic.html">Jan Marjanovic</a>
        </li>
        <li class="published" title="2014-11-06T22:00:00+01:00">
          on&nbsp;Thu 06 November 2014
        </li>

	</ul>
<p>Category: <a href="../tag/theremin.html">Theremin</a>, </p>
</div><!-- /.post-info -->          <div class="highlight"><pre><span></span><code><span class="c1"># -*- coding: utf-8 -*-</span>
<span class="sd">&quot;&quot;&quot;</span>
<span class="sd">Created on Sun Nov  2 14:54:35 2014</span>

<span class="sd">@author: jan</span>
<span class="sd">&quot;&quot;&quot;</span>

<span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
<span class="kn">from</span> <span class="nn">mpl_toolkits.mplot3d</span> <span class="kn">import</span> <span class="n">Axes3D</span>
<span class="kn">from</span> <span class="nn">matplotlib</span> <span class="kn">import</span> <span class="n">cm</span>
<span class="kn">import</span> <span class="nn">matplotlib.pyplot</span> <span class="k">as</span> <span class="nn">plt</span>
<span class="kn">import</span> <span class="nn">scipy.constants</span> <span class="k">as</span> <span class="nn">const</span>

<span class="c1">#   b, y</span>
<span class="c1"># A</span>
<span class="c1"># |</span>
<span class="c1"># +----&gt; a, x</span>
<span class="n">a</span> <span class="o">=</span> <span class="mf">0.3</span>         <span class="c1"># x dimension of space</span>
<span class="n">b</span> <span class="o">=</span> <span class="mf">1.0</span>         <span class="c1"># y dimension of space</span>

<span class="n">Na</span> <span class="o">=</span> <span class="mi">50</span>         <span class="c1"># number of nodes in x</span>
<span class="n">Nb</span> <span class="o">=</span> <span class="mi">100</span>        <span class="c1"># number of nodes in y</span>

<span class="n">x_ant</span> <span class="o">=</span> <span class="mf">0.1</span>     <span class="c1"># antenna x position</span>
<span class="n">y_ant</span> <span class="o">=</span> <span class="mf">0.5</span>     <span class="c1"># antenna y position</span>
<span class="n">l_ant</span> <span class="o">=</span> <span class="mf">0.5</span>     <span class="c1"># antenna length</span>
<span class="n">V_ant</span> <span class="o">=</span> <span class="mf">1.0</span>     <span class="c1"># antenna voltage (capacitance is Q/V, so we set V to 1)</span>

<span class="n">x_man</span> <span class="o">=</span> <span class="mf">0.2</span>     <span class="c1"># man x position</span>
<span class="n">y_man</span> <span class="o">=</span> <span class="mf">0.4</span>     <span class="c1"># man y position</span>
<span class="n">l_man</span> <span class="o">=</span> <span class="mf">0.8</span>     <span class="c1"># man length (height)</span>
<span class="n">V_man</span> <span class="o">=</span> <span class="mf">0.0</span>     <span class="c1"># man voltage</span>

<span class="n">l_hand</span> <span class="o">=</span> <span class="mf">0.06</span>   <span class="c1"># hand length   &lt;--- VERY IMPORTANT PARAMETER</span>
<span class="n">x_hand</span> <span class="o">=</span> <span class="n">x_man</span><span class="o">-</span><span class="n">l_hand</span><span class="o">/</span><span class="mf">2.0</span>
<span class="n">y_hand</span> <span class="o">=</span> <span class="mf">0.6</span>    <span class="c1"># hand y position</span>
<span class="n">V_hand</span> <span class="o">=</span> <span class="n">V_man</span>  <span class="c1"># hand voltage has same potential as man</span>


<span class="n">_err_ppm</span> <span class="o">=</span> <span class="mi">1000</span>   <span class="c1"># max diff from two consecutive iterations to stop the loop</span>

<span class="c1">###############################################################################</span>
<span class="c1"># Setup</span>
<span class="c1">#   We create 2D grid of nodes with their voltage and save if they are set</span>
<span class="c1">#   to fixed voltage (e.g. grounded) or is there a field which should be </span>
<span class="c1">#   caculated</span>

<span class="n">V</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">Na</span><span class="p">,</span><span class="n">Nb</span><span class="p">))</span>
<span class="n">fixed</span> <span class="o">=</span> <span class="p">[[</span><span class="kc">False</span><span class="p">]</span><span class="o">*</span><span class="n">Nb</span> <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Na</span><span class="p">)]</span>

<span class="k">def</span> <span class="nf">placeObject</span><span class="p">(</span><span class="n">x</span><span class="p">,</span><span class="n">y</span><span class="p">,</span><span class="n">l</span><span class="p">,</span><span class="n">vert</span><span class="p">,</span><span class="n">V0</span><span class="p">):</span>
    <span class="k">global</span> <span class="n">V</span>
    <span class="k">if</span> <span class="n">vert</span><span class="p">:</span> 
    <span class="n">Nx</span>       <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">x</span><span class="o">/</span><span class="n">a</span><span class="o">*</span><span class="n">Na</span><span class="p">)</span>
    <span class="n">Ny_start</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">y</span><span class="o">-</span><span class="n">l</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">/</span><span class="n">b</span><span class="o">*</span><span class="n">Nb</span><span class="p">)</span> 
    <span class="n">Ny_stop</span>  <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">y</span><span class="o">+</span><span class="n">l</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">/</span><span class="n">b</span><span class="o">*</span><span class="n">Nb</span><span class="p">)</span>

    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Ny_start</span><span class="p">,</span> <span class="n">Ny_stop</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
        <span class="n">V</span><span class="p">[</span><span class="n">Nx</span><span class="p">][</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="n">V0</span>
        <span class="n">fixed</span><span class="p">[</span><span class="n">Nx</span><span class="p">][</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>
    <span class="k">else</span><span class="p">:</span>
    <span class="n">Nx_start</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">x</span><span class="o">-</span><span class="n">l</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">/</span><span class="n">a</span><span class="o">*</span><span class="n">Na</span><span class="p">)</span>
    <span class="n">Nx_stop</span>  <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">x</span><span class="o">+</span><span class="n">l</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">/</span><span class="n">a</span><span class="o">*</span><span class="n">Na</span><span class="p">)</span>
    <span class="n">Ny</span>       <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">y</span><span class="o">/</span><span class="n">b</span><span class="o">*</span><span class="n">Nb</span><span class="p">)</span>

    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Nx_start</span><span class="p">,</span> <span class="n">Nx_stop</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
        <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">Ny</span><span class="p">]</span> <span class="o">=</span> <span class="n">V0</span>            
        <span class="n">fixed</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">Ny</span><span class="p">]</span> <span class="o">=</span> <span class="kc">True</span>


<span class="n">placeObject</span><span class="p">(</span><span class="n">x_ant</span><span class="p">,</span>  <span class="n">y_ant</span><span class="p">,</span>  <span class="n">l_ant</span><span class="p">,</span>  <span class="kc">True</span><span class="p">,</span>  <span class="n">V_ant</span><span class="p">)</span>
<span class="n">placeObject</span><span class="p">(</span><span class="n">x_man</span><span class="p">,</span>  <span class="n">y_man</span><span class="p">,</span>  <span class="n">l_man</span><span class="p">,</span>  <span class="kc">True</span><span class="p">,</span>  <span class="n">V_man</span><span class="p">)</span>
<span class="n">placeObject</span><span class="p">(</span><span class="n">x_hand</span><span class="p">,</span> <span class="n">y_hand</span><span class="p">,</span> <span class="n">l_hand</span><span class="p">,</span> <span class="kc">False</span><span class="p">,</span> <span class="n">V_hand</span><span class="p">)</span>

<span class="c1">###############################################################################</span>
<span class="c1"># Calculation</span>
<span class="c1">#   We are solving 1st Maxwell equation:</span>
<span class="c1">#      div D = ro -&gt; div E = ro/eps --&gt; div E = 0 </span>
<span class="c1">#</span>
<span class="c1">#   Combined with definiton of voltage potential:</span>
<span class="c1">#      E = - grad V</span>
<span class="c1">#</span>
<span class="c1">#   Results in Laplace equation:</span>
<span class="c1">#      div grad V = 0</span>


<span class="n">m</span><span class="p">,</span> <span class="n">n</span> <span class="o">=</span> <span class="n">V</span><span class="o">.</span><span class="n">shape</span>

<span class="k">def</span> <span class="nf">iterate</span><span class="p">():</span>
    <span class="k">global</span> <span class="n">V</span>
    <span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">m</span><span class="p">):</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
        <span class="k">if</span> <span class="n">i</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">and</span> <span class="n">j</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="mi">0</span><span class="p">])</span> <span class="o">/</span> <span class="mf">2.0</span>
        <span class="k">elif</span> <span class="n">i</span> <span class="o">==</span> <span class="n">m</span><span class="o">-</span><span class="mi">1</span> <span class="ow">and</span> <span class="n">j</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">1</span><span class="p">])</span> <span class="o">/</span> <span class="mf">2.0</span>
        <span class="k">elif</span> <span class="n">i</span> <span class="o">==</span> <span class="mi">0</span> <span class="ow">and</span> <span class="n">j</span> <span class="o">==</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">])</span> <span class="o">/</span> <span class="mf">2.0</span>        
        <span class="k">elif</span> <span class="n">i</span> <span class="o">==</span> <span class="n">m</span><span class="o">-</span><span class="mi">1</span> <span class="ow">and</span> <span class="n">j</span> <span class="o">==</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">])</span> <span class="o">/</span> <span class="mf">2.0</span>   
        <span class="k">elif</span> <span class="n">i</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="mi">0</span><span class="p">][</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">])</span> <span class="o">/</span> <span class="mf">3.0</span> 
        <span class="k">elif</span> <span class="n">i</span> <span class="o">==</span> <span class="n">m</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">])</span> <span class="o">/</span> <span class="mf">3.0</span>      
        <span class="k">elif</span> <span class="n">j</span> <span class="o">==</span> <span class="mi">0</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">][</span><span class="mi">0</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="mi">1</span><span class="p">])</span> <span class="o">/</span> <span class="mf">3.0</span> 
        <span class="k">elif</span> <span class="n">j</span> <span class="o">==</span> <span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">])</span> <span class="o">/</span> <span class="mf">3.0</span>    
        <span class="k">else</span><span class="p">:</span>
            <span class="n">V_T</span> <span class="o">=</span> <span class="p">(</span><span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">-</span><span class="mi">1</span><span class="p">]</span> <span class="o">+</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">])</span> <span class="o">/</span> <span class="mf">4.0</span>

        <span class="k">if</span> <span class="ow">not</span> <span class="n">fixed</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]:</span>
            <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">V_T</span>


<span class="n">e</span> <span class="o">=</span> <span class="p">[]</span>
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="mi">1000</span><span class="p">):</span>
    <span class="n">Vprev</span> <span class="o">=</span> <span class="n">V</span><span class="o">.</span><span class="n">copy</span><span class="p">()</span>
    <span class="n">iterate</span><span class="p">()</span>
    <span class="n">ei</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sum</span><span class="p">(</span><span class="n">Vprev</span> <span class="o">-</span> <span class="n">V</span><span class="p">)</span>
    <span class="n">e</span><span class="o">.</span><span class="n">append</span><span class="p">(</span><span class="n">ei</span><span class="p">)</span>
    <span class="k">if</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">abs</span><span class="p">(</span><span class="n">ei</span><span class="o">/</span><span class="n">Na</span><span class="o">/</span><span class="n">Nb</span><span class="p">)</span> <span class="o">&lt;</span> <span class="n">_err_ppm</span><span class="o">*</span><span class="mf">1e-6</span><span class="p">):</span>
    <span class="nb">print</span> <span class="s2">&quot;Done&quot;</span>
    <span class="k">break</span>


<span class="c1">###############################################################################</span>
<span class="c1"># Calc C</span>
<span class="c1">#   Electric field is calculated as E = - grad V</span>
<span class="c1">#   Charge is calculated from Gauss equation in integral form</span>

<span class="n">Ex</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">Na</span><span class="p">,</span><span class="n">Nb</span><span class="p">))</span>
<span class="n">Ey</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">Na</span><span class="p">,</span><span class="n">Nb</span><span class="p">))</span>
<span class="n">Eabs</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">zeros</span><span class="p">((</span><span class="n">Na</span><span class="p">,</span><span class="n">Nb</span><span class="p">))</span> 
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">m</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
    <span class="k">for</span> <span class="n">j</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="o">-</span><span class="mi">1</span><span class="p">):</span>
    <span class="n">Ex</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="o">+</span><span class="mi">1</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">-</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">a</span><span class="o">/</span><span class="n">Na</span><span class="p">)</span>
    <span class="n">Ey</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="p">(</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">-</span> <span class="n">V</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="p">)</span> <span class="o">/</span> <span class="p">(</span><span class="n">b</span><span class="o">/</span><span class="n">Nb</span><span class="p">)</span>
    <span class="n">Eabs</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="n">Ex</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">Ex</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="o">+</span><span class="n">Ey</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">]</span><span class="o">*</span><span class="n">Ey</span><span class="p">[</span><span class="n">i</span><span class="p">][</span><span class="n">j</span><span class="p">])</span>



<span class="n">Nx</span>       <span class="o">=</span> <span class="nb">int</span><span class="p">(</span><span class="n">x_ant</span><span class="o">/</span><span class="n">a</span><span class="o">*</span><span class="n">Na</span><span class="p">)</span>
<span class="n">Ny_start</span> <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">y_ant</span><span class="o">-</span><span class="n">l_ant</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">/</span><span class="n">b</span><span class="o">*</span><span class="n">Nb</span><span class="p">)</span> 
<span class="n">Ny_stop</span>  <span class="o">=</span> <span class="nb">int</span><span class="p">((</span><span class="n">y_ant</span><span class="o">+</span><span class="n">l_ant</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span><span class="o">/</span><span class="n">b</span><span class="o">*</span><span class="n">Nb</span><span class="p">)</span>


<span class="c1"># Gauss </span>
<span class="n">Q</span> <span class="o">=</span> <span class="mi">0</span>

<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">Ny_start</span><span class="p">,</span> <span class="n">Ny_stop</span><span class="o">+</span><span class="mi">1</span><span class="p">):</span>
    <span class="n">Q</span> <span class="o">+=</span> <span class="n">Eabs</span><span class="p">[</span><span class="n">Nx</span><span class="o">+</span><span class="mi">1</span><span class="p">][</span><span class="n">i</span><span class="p">]</span>
    <span class="n">Q</span> <span class="o">+=</span> <span class="n">Eabs</span><span class="p">[</span><span class="n">Nx</span><span class="o">-</span><span class="mi">1</span><span class="p">][</span><span class="n">i</span><span class="p">]</span>

<span class="n">C</span> <span class="o">=</span> <span class="n">const</span><span class="o">.</span><span class="n">epsilon_0</span> <span class="o">*</span> <span class="n">Q</span>

<span class="c1">###############################################################################</span>
<span class="c1"># Plot</span>
<span class="n">X</span><span class="p">,</span><span class="n">Y</span> <span class="o">=</span> <span class="n">np</span><span class="o">.</span><span class="n">meshgrid</span><span class="p">(</span><span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">b</span><span class="p">,</span> <span class="n">Nb</span><span class="p">),</span> <span class="n">np</span><span class="o">.</span><span class="n">linspace</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="n">a</span><span class="p">,</span> <span class="n">Na</span><span class="p">))</span>
<span class="n">fig</span> <span class="o">=</span> <span class="n">plt</span><span class="o">.</span><span class="n">figure</span><span class="p">(</span><span class="n">figsize</span><span class="o">=</span><span class="p">(</span><span class="mi">8</span><span class="p">,</span><span class="mi">6</span><span class="p">))</span>
<span class="n">ax</span> <span class="o">=</span> <span class="n">fig</span><span class="o">.</span><span class="n">gca</span><span class="p">(</span><span class="n">projection</span><span class="o">=</span><span class="s1">&#39;3d&#39;</span><span class="p">)</span>

<span class="n">ax</span><span class="o">.</span><span class="n">view_init</span><span class="p">(</span><span class="o">-</span><span class="mi">90</span><span class="p">,</span><span class="mi">0</span><span class="p">)</span>
<span class="n">ax</span><span class="o">.</span><span class="n">plot_surface</span><span class="p">(</span><span class="n">X</span><span class="p">,</span> <span class="n">Y</span><span class="p">,</span> <span class="n">V</span><span class="p">,</span> <span class="n">rstride</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">cstride</span><span class="o">=</span><span class="mi">1</span><span class="p">,</span> <span class="n">cmap</span><span class="o">=</span><span class="n">cm</span><span class="o">.</span><span class="n">rainbow</span><span class="p">,</span>
    <span class="n">linewidth</span><span class="o">=</span><span class="mi">0</span><span class="p">,</span> <span class="n">antialiased</span><span class="o">=</span><span class="kc">False</span><span class="p">)</span>

<span class="n">ax</span><span class="o">.</span><span class="n">text</span><span class="p">(</span><span class="mf">0.95</span><span class="p">,</span> <span class="mf">2.5</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="s1">&#39;C = &#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">(</span><span class="nb">round</span><span class="p">(</span><span class="n">C</span><span class="o">*</span><span class="mf">1e9</span><span class="p">,</span><span class="mi">4</span><span class="p">))</span> <span class="o">+</span> <span class="s2">&quot; nF&quot;</span><span class="p">,</span>
    <span class="n">verticalalignment</span><span class="o">=</span><span class="s1">&#39;bottom&#39;</span><span class="p">,</span> <span class="n">horizontalalignment</span><span class="o">=</span><span class="s1">&#39;right&#39;</span><span class="p">,</span>
    <span class="n">transform</span><span class="o">=</span><span class="n">ax</span><span class="o">.</span><span class="n">transAxes</span><span class="p">,</span> <span class="n">fontsize</span><span class="o">=</span><span class="mi">12</span><span class="p">)</span>

<span class="n">ax</span><span class="o">.</span><span class="n">set_title</span><span class="p">(</span><span class="s1">&#39;Distance = &#39;</span> <span class="o">+</span> <span class="nb">str</span><span class="p">((</span><span class="n">x_hand</span><span class="o">-</span><span class="n">l_hand</span><span class="o">/</span><span class="mf">2.0</span><span class="p">)</span> <span class="o">-</span> <span class="n">x_ant</span> <span class="p">)</span> <span class="o">+</span> <span class="s1">&#39; cm&#39;</span><span class="p">)</span>
</code></pre></div>
        </div><!-- /.entry-content -->
<a href="https://twitter.com/share" class="twitter-share-button" data-count="horizontal" data-via="janmarjanovic">Tweet</a><script type="text/javascript" src="https://platform.twitter.com/widgets.js"></script><br/><br/>

      </article>
    </div>
</section>
        <section id="extras" >
        
        </section><!-- /#extras -->
	
        <footer id="contentinfo" >
                <address id="about" class="vcard ">
                Proudly powered by <a href="https://getpelican.com/" target="_blank">Pelican</a>, which takes
                great advantage of <a href="https://python.org" target="_blank">Python</a>.
		
                </address><!-- /#about -->
		

                
        </footer><!-- /#contentinfo -->

<script>
  (function(i,s,o,g,r,a,m){i['GoogleAnalyticsObject']=r;i[r]=i[r]||function(){
  (i[r].q=i[r].q||[]).push(arguments)},i[r].l=1*new Date();a=s.createElement(o),
  m=s.getElementsByTagName(o)[0];a.async=1;a.src=g;m.parentNode.insertBefore(a,m)
  })(window,document,'script','//www.google-analytics.com/analytics.js','ga');

  ga('create', 'UA-56555055-1', 'auto');
  ga('send', 'pageview');

</script></body>
</html>